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Abstract 

On the basis of the extensive air shower (EAS) data obtained by the GAMMA 
experiment, the energy spectra and elemental composition of the primary cosmic 
rays are derived in the 10^ — 10^ TeV energy range. The reconstruction of the primary 
energy spectra is carried out using an EAS inverse approach in the framework of 
the SIBYLL2.1 and QGSJETOl interaction models and the hypothesis of power- 
law primary energy spectra with rigidity-dependent knees. The energy spectra of 
primary H, He, O-like and Fe-like nuclei obtained with the SIBYLL interaction 
model agree with corresponding extrapolations of the balloon and satellite data to 
~ 10^ TeV energies. The energy spectra obtained from the QGSJET model show a 
predominantly proton composition in the knee region. The rigidity-dependent knee 
feature of the primary energy spectra for each interaction model is displayed at the 
following rigidities: Er ~ 2500 ± 200 TV (SIBYLL) and Er ~ 3100 - 4200 TV 
(QGSJET). 

All the results presented are derived taking into account the detector response, 
the reconstruction uncertainties of the EAS parameters, and fluctuations in the EAS 
development. 

Key words: Cosmic rays, energy spectra, composition, extensive air showers 
PACS: 96.40.Pq, 96.40.De, 96.40.-z, 98.70.Sa 



* Corresponding author: 

Email address: samvelSyerphi . am (S.V. Ter- Antonyan) . 



Preprint submitted to Elsevier Science 



5 February 2008 



1 Introduction 



The investigation of the energy spectra and elemental composition of primary 
cosmic rays in the knee region (10^ — 10^ TeV) remains one of the intriguing 
problems of modern high energy cosmic-ray physics. Despite the fact that 
these investigations have been carried out for more than half a century, the 
data on the elemental primary energy spectra at energies E > 10^ TeV still 
need improvement. The high statistical accuracies of recent EAS experiments 
[1,2,3,4] have confirmed the presence of a bend in the all-particle primary 
energy spectrum at around 3 • 10^ TeV (called the "knee") from an overall 
spectrum oc E~'^-'^ below the knee to oc E~^-^ beyond the knee, and a change in 
composition toward heavier species with increasing energy in the 10^ — 10^ TeV 
region. However, separating the primary energy spectra of elemental groups 
remains difficult, both due to uncertainties in the interaction model and the 
uncertainties associated with the solutions to the EAS inverse problem [5,6]. 

One of the most studied class of models for the origin of cosmic rays in this 
energy region, which assumes that supernova remnants are their main source, 
predicts rigidity-dependent primary energy spectra in the knee region ([7,8] 
and references therein) . Other astrophysical models for the origin of the knee, 
such as Galactic propagation effects [9,10] also predict rigidity-dependent spec- 
tra. Such energy spectra of primary nuclei with rigidity-dependent knees can 
approximately describe the observed EAS size spectra in the 10^ — 10^ TeV en- 
ergy region in the framework of conventional interaction models [11,12,13,14,15] 
However, an alternative class of models predicts mass- dependent knees (see 
[16] and references therein for a recent review of models of the origin of the 
knee). In the present analysis, we will assume a rigidity-dependent knee; the 
appropriateness of this hypothesis will be briefly examined in the discussion 
section. 

The GAMMA facility (Fig. 1) was designed at the beginning of the 1990's 
in the framework of the ANI experiment [17] and the first results of EAS 
investigations were presented in [18,19,20,21]. The main characteristic features 
of the GAMMA experiment are the mountain location, the symmetric location 
of the EAS detectors, and the underground muon scintillation carpet which 
detects the EAS muon component with energy E^ > 5 GeV. 

Here, a description of the GAMMA facihty, the main results of investigation 

during 2002-2004 [20,21,22] and evaluations of primary energy spectra in the 
knee region are presented in comparison with the corresponding simulated data 
in the framework of the SIBYLL [23] and QGSJET [24] interaction models. 
Prehminary results have already been presented in [25,26,27]. 
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Fig. 1. Diagrammatic layout of the GAMMA facility. 



2 GAMMA experiment 

The GAMMA installation [18,19,20,21,22] is a ground-based array of 33 sur- 
face particle detection stations and 150 underground muon detectors, located 
on the south side of Mount Aragats in Armenia. The elevation of the GAMMA 
facility is 3200 m above sea level, which corresponds to 700 g/cm^ of atmo- 
spheric depth. A diagrammatic layout is shown in Fig. 1. 
The surface stations of the EAS array are located on 5 concentric circles of 
radii ~20, 28, 50, 70 and 100 m, and each station contains 3 square plastic scin- 
tillation detectors with the following dimensions: 1 x 1 x 0.05 m^. Each of the 
central 9 stations contains an additional (4th) small scintillator with dimen- 
sions 0.3 X 0.3 X 0.05 m'^ (Fig. 1) for high particle density (3> 10^ particles/m^) 
measurements. 

A photomultiplier tube is positioned on the top of the aluminum casing cover- 
ing each scintillator. One of the three detectors of each station is examined by 
two photomultipliers, one of which is designed for fast timing measurements. 
150 underground muon detectors (muon carpet) are compactly arranged in 
the underground hall under 2.3 kg/cm^ of concrete and rock. The scintillator 
dimensions, casings and photomultipliers are the same as in the EAS surface 
detectors. 



2.1 Detector system and triggering 

The output voltage of each photomultiplier is converted into a pulse burst 
by a logarithmic ADC and transmitted to a CAMAC array, where the corre- 
sponding electronic counters produce a digital number ("code") of pulses in 
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the burst. Four inner ("trigger") stations at a radius of 20m are monitored 
by a coincidence circuit. If at least two scintillators of each trigger station 
each detect more than 3 particles, the information from all detectors are then 
recorded along with the time between the master trigger pulse and the pulses 
from all fast-timing detectors. The given trigger condition provides EAS de- 
tections with an EAS size threshold Nch > (0.5 — 1) • 10^ for a location of the 
EAS core within the i? < 50 m circle. The shower size thresholds for 100% 
shower detection efficiency are equal to Nch = 3 • 10^ and N^h = 5 • 10^ for 
EAS core locations within i? < 25 m and i? < 50 m respectively [18]. 
Before being placed on the scintillator casing, all photomultiphers were tested 
by a test bench using a luminodiode method where the corresponding param- 
eters of the logarithmic ADC and the upper limit ((0.5 — 1) ■ 10^) [28] of the 
particle density measurement ranges were determined for each detector. The 
number of charged particles (rij) passing through the i-th scintillator is calcu- 
lated using a logarithmic transformation: Inrij — {C — Co)/d [28], where the 
scale parameter ~ (9 — 10) ± 0.35 is determined for each detector by the 
test bench, < C < 2^ — 1 is the output digital code from the CAMAC array 
corresponding to the energy deposit of n charged particles into the scintilla- 
tor, and Co — (5 — 6) ± 0.25 is equal to the mode of the background single 
particle digital code spectra (Section 2.4). The time delay is estimated by the 
pair-delay method [29] to give a time resolution of about 4 — 5 ns. 

2.2 Reconstruction of EAS parameters 

The EAS zenith angle (6) is estimated on the basis of the shower front arrival 
times measured by the 33 fast-timing surface detectors, applying a maximum 
likelihood method and the flat-front approach [29,30]. The corresponding un- 
certainty was tested by Monte Carlo simulations and is equal to a{9) ~ 1.5° 
[18]. The reconstruction of the EAS size {N^h), shower age (s) and core coor- 
dinates {xo,yo) is performed based on the Nishimura-Kamata-Creisen (NKG) 
approximation to the measured charged particle densities {{ni}, i = 1, . . . ,m), 
using minimization to estimate Xo,yo and a maximum likelihood method 
to estimate N^h, taking into account the measurement errors. Gamma-quanta 
conversions in the scintillator and housing were taken into account in the es- 
timates of Nch (Section 2.3). 

The logarithmic transformation L{ni) = Inn^ — (1/m) 5]]lnnj for Ui ^ en- 
ables an analytical solution for the EAS age parameter (s) using mini- 
mization [30,31]. Unbiased (< 5%) estimations of Nch,s,Xo and i/q shower 
parameters are obtained for Nch > 5 • 10^, 0.3 < s < 1.6, 9 < 30° and dis- 
tances i? < 25 m from the shower core to the center of the EAS array. The 
shower age parameter (s) is estimated from the surface scintillators located 
inside a 7 m < i?j < 80 m ring area around the shower core (Section 2.3). 
The EAS detection efficiency (P^) and corresponding accuracies are derived 
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from mimic shower simulations taking into account the EAS fluctuations and 
measurement errors (Section 2.4) and are equal to: = 100%, ANf.fjNch — 
0.1, As ~ 0.05, Ax and Ay ~ 0.5 — 1 m. These results were also checked with 
CORSIKA [32] simulated EAS (Section 2.3) and depend slightly on shower 
core location for i? < 50 m. 

The reconstruction of the total number of EAS muons (Nfj) from the de- 
tected muon densities {{n^j},j = 1, • • • , 150) in the underground muon hall 
is carried out by restricting the distance to < 50 m from the shower core 
(the so-called EAS "truncated" muon size [18,33]) and using the approxima- 
tion to the muon lateral distribution function [18,34]: Pij,{r) — cNn{Rfj_ < 
50m) exp (— r/ro)/(r/ro)°'''', where tq = 80 m [35] and c = l/27r J^^ p{r)rdr. 
The EAS truncated muon size N^{R^ < 50m) is estimated at known (from 
the EAS surface array) shower core coordinates in the underground muon hall. 
Unbiased estimations for muon size are obtained for > 10^ using a max- 
imum likelihood method and assuming Poisson fluctuations in the detected 
muon numbers. The reconstruction accuracies of the truncated muon shower 
sizes are equal to AN^/N^ ^ 0.2 — 0.35 at ~ 10^ — 10^ respectively. 
It should be noted that the detected muons in the underground hall are always 
accompanied by the electron-positron equilibrium spectrum which is produced 
when muons pass through the matter (2300 g/cm^) over the scintillation car- 
pet; this is taken into account in our results (Section 3.2). 

2.3 Detector response 

The GAMMA detector response taking into account the EAS 7-quanta con- 
tribution was computed by EAS simulations using the CORSIKA 6.031 code 
[32] (NKG and EGS modes, GHEISHA2002) with the QGSJETOl [24] and 
SIBYLL 2.1 [23] interaction models for 4 types {A = H,He,0,Fe) of pri- 
mary nuclei. Each EAS particle (7, e, p, h) obtained from CORSIKA (EGS 
mode) at the observation level was examined by passing through the steel 
casing (1.5 mm) of the detector station and then through the corresponding 
scintillator. The pair production and Compton scattering processes were addi- 
tionally simulated in the case of EAS 7-quanta. The resulting energy deposit 
in the scintillator was converted to an ADC code and inverse-decoded into a 
number of "detected" charged particles taking into account all uncertainties 
of the ADC parameters (Co, d) and fluctuations in the light collected by the 
photomultipliers (cx; ~ 0.25/v^). 

Using the simulation scenario above, 200 EAS events with shower size thresh- 
old Nch > 5 • 10^ were simulated with CORSIKA simultaneously in the EGS 

and NKG modes for each of the A = He, O and Fe primary nuclei, with 
logarithm-uniform energy spectra in the 10^ — 10^ TeV energy range. The 
computation of the charged particle densities in surface detectors in the NKG 
mode was performed by applying two-dimensional interpolations of the corre- 
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spending shower electron (and positron) density matrix from CORSIKA [32], 
along with the individual EAS muons and hadrons. 

A ~ 5% agreement between the EGS (including the EAS 7-quanta contribu- 
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Fig. 2. Biases (a) and standard deviations (b) versus EAS size, and distributions 

of biases in shower size (c) and shower age parameter (d) for the SIBYLL (circle 
symbols, lines) and QGSJET (square symbols) interaction models for primary H 
(empty symbols, solid lines) and Fe (filled symbols, dashed lines) nuclei. 



tion) and NKG simulated EAS data was attained for sm E^'^l MeV kinetic 
energy threshold of shower electrons (and positrons) in the NKG mode, con- 
sidering only the 7m < Ri < 80m ring area used in the determination of the 
shower age parameter. Thus the underestimation of the EAS particle density 
due to the threshold of the detected energy deposit {Ed ~ 8 MeV [18,25]) in 
the scintillators is compensated by the EAS 7-quanta contribution. 
The corresponding biases 

. ,_ N,h{Ee = lMey,NKG) 



and standard deviations (T(5jVcft) versus the reconstructed EAS size (Nch) are 
shown in Fig. 2 (a) and (b) respectively, for the SIBYLL (circle symbols) 
and QGSJET (square symbols) interaction models and for primary H (empty 
symbols) and Fe nuclei (filled symbols). The distributions of the biases in 
reconstructed EAS sizes (^at^^) and shower age parameters 

Ss{A) = s{Ee = 1 MeV, NKG) - s{Ed, 7, EGS) 
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arc shown in Fig. 2 (c) and (d) respectively, for a shower size threshold 
Nch > 5-10^, the SIBYLL interaction model, and primary H (solid lines) 
and Fe (dashed lines) nuclei. 

The observed (~ 5%) biases in 5^^^ (Fig. 2a) for the 4 kinds of primary nu- 
clei depend only weakly on the interaction model (< 5%) and zenith angles 

(< 3% for 6 < 30°), and the biases in age parameter dg can be considered neg- 
ligible. The NKG-mode simulated sizes were divided by the estimated biases 
Sn^^{A, Nch) in the reconstruction of the primary energy spectra (Section 3.1). 



2.4 Measurement errors and density spectra 

The close disposition of the k = 1, 2, 3 scintillators in each of the (i-th) detector 
stations of the GAMMA surface array enables a calibration of the measure- 
ment error using the detected EAS data. The measured and simulated particle 
density discrepancies (n^ — p)/p versus the average value p = (l/3)X^nfe for 
distances i?j > 10 m from the shower core are shown in Fig. 3 (circle symbols), 
and are completely determined by Poisson fluctuations (at i?j 3> 1 m ) and 
the measurement errors. The agreement between the measured and simulated 
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Fig. 3. Particle density discrepancies (circle symbols) and measurement error for a 
single detector (square symbols) versus charged-particle density. 



dependences enables the extraction of the actual measurement errors of the 
GAMMA detectors. The corresponding results, obtained from the simulations 
without Poisson fluctuations, are shown in Fig. 3 (square symbols). 
The background omni-directional single particle spectra (in units of ADC 
code) detected by GAMMA surface scintillators in 78 s of operation time 



Q. 



O 1 



0) 

E 



0.5 



7 



arc shown in Fig. 4 (dotted lines). The background single particle spectra de- 
tected by underground muon scintillators have the same shape but about 10 
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Fig. 4. Background single particle spectra of 15 surface detectors (dotted lines). 
The symbols (solid line) are the expected spectra taking into account (without) 
measurement errors. 



times lower intensities. These spectra (pulse height distributions) along with 
the known zenith angle distributions and composition (~ 40%e, 50%/i) of the 
background charged particles at the observation level [36] are used for the 
operative determination of the ADC parameters (Co) for each experimental 
run. The symbols and solid lines in Fig. 4 display the corresponding expected 
spectra obtained by CORSIKA (EGS) simulation, without errors (solid line) 
and taking into account the measurement errors (symbols) respectively. The 
minimal primary energy in the simulation of the background particle spectra 
was determined by the 7.6 GV geomagnetic rigidity cutoff in Armenia. 
Because the effective primary energies responsible for the single particle spec- 
tra at the observation level of 700 g/cw? are around 100 GeV, and this energy 
range has been studied by direct measurements in balloon and satellite exper- 
iments, the primary energy spectra and elemental composition in the Monte 
Carlo simulation were taken from power-law approximations to the direct mea- 
surement data [37]. It should be noted that the expected single particle spectra 
at these energies arc practically the same for the QGSJET and SIBYLL in- 
teraction models, because most of the interactions occur in the energy range 
where accelerator data are used. 

Fig. 5 (symbols) displays the EAS charged particle density spectra measured 
by the surface detectors (left panel) and underground muon detectors (right 
panel) at i?j < 50 m with different EAS size thresholds: N^h > 5-10''', N^h > 10^ 
(and additionally Nch > 2 • 10® for the muon density spectra). The showers 
were selected with 9 < 30° and shower core location in the i? < 25 m range 
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Fig. 5. Detected (symbols) and expected (lines) particle density spectra measured 
by the surface scintillators (left panel) and underground muon scintillators (right 
panel). 



from the center of the GAMMA facility (Fig. 1). The corresponding expected 
spectra (lines) for different interaction models are also shown in Fig. 5. The 
primary energy spectra and elemental composition for these simulations were 
those obtained in the combined approximation solution to the EAS inverse 
problem (Section 3.3). There is good agreement between the expected and 
observed data for the surface array over the full measurement range (about 
four orders of magnitude). However, agreement of the detected muon den- 
sity spectra with the expected ones is attained only in the Nch < 10^ range. 
The observed discrepancies for the muon density spectra at Nch > 10^ are 
unaccounted for at present, and will require subsequent investigations. 



2.5 EAS data set 



The data set analyzed in this paper was obtained over 6.19 • 10^ s of operating 
live time of the GAMMA facility, from 2002 to 2004. Showers were selected 
for analysis with the following criteria: Nch > 5 ■ 10^ , R < 25 m, 6 < 30°, 
0.3 < s < 1.6, X^^Nch) /m < 3 and 'X^{N^)/m < 3 (where m is the number of 
scintillators with non-zero signal), yielding a total data set of 1.9 • 10^ selected 
showers. The selected measurement range provided 100% EAS detection effi- 
ciency (Section 2.2) and similar conditions for the reconstruction of the shower 
lateral distribution functions. 

The measured variable distributions used in the combined approximation ap- 
proach to the EAS inverse problem (Section 3.3) are shown in Figs. 6-11 (sym- 
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EAS size, N 



Fig. 6. EAS size spectra for 3 zenith angle intervals (symbols) and corresponding 
expected spectra according to the SIBYLL (solid lines) and QGSJET (dashed lines) 
interaction models. 




Fig. 7. Normalized EAS truncated muon size spectra for 3 zenith angle intervals 
(symbols) and corresponding expected spectra for the SIBYLL (solid lines) and 
QGSJET (dashed lines) interaction models. 



bols). All lines and shaded areas in these figures correspond to the expected 
spectra computed on the basis of the EAS inverse problem solution in the 
framework of the SIBYLL and QGSJET interaction models. These expected 
(forward folded) spectra are computed by Monte-Carlo integration (Section 
3.1) using the simulated EAS database, which results in the statistical fiuctu- 
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Fig. 8. EAS size spectra (symbols) for different truncated muon size thresholds 
(Nfj^) and corresponding expected spectra according to the SIBYLL (solid lines) 
and QGSJET (dashed lines) interaction models. 




N^(R^<50m) 



Fig. 9. Normalized EAS truncated muon size spectra (symbols) for different shower 
size thresholds (Nch) and corresponding expected spectra according to the SIBYLL 
(solid lines) and QGSJET (dashed lines) interaction models. 



ations evident in many of these predicted spectra. 

The EAS size spectra {N^h ' d-F(6') / dNch) for three zenith angle intervals 
are shown in Fig. 6. The EAS truncated muon size spectra in the same zenith 
angle intervals are shown in Fig. 7; these spectra are normalized per unit 
shower with N^h > 5 ■ 10^ and ^ < 30°. The EAS size spectra for 9 < 30° 
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Fig. 10. Dependence of the average EAS age parameter on EAS size (symbols) along 
with expected values for the SIBYLL (solid line) and QGSJET (dashed line) inter- 
action models. The dotted and dash-dotted lines correspond to expected values for 
primary hydrogen and iron nuclei, respectively, for the SIBYLL interaction model. 




Fig. 11. Average EAS truncated muon size A^^ versus EAS size N^h (symbols). The 
lines and shaded area are the expected dependences for the SIBYLL and QGSJET 
interaction models and primary p, Fe and mixed (Table 1) compositions. 



and different thresholds in EAS truncated muon size are shown in Fig. 8. The 
normalized EAS truncated muon size spectra for different EAS size thresholds 
are shown in Fig. 9. Fig. 10 displays the dependence of the average EAS age 
parameter on EAS size s{Nch)- Fig. 11 shows the observed N^{Nch) depen- 
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dence and the corresponding expected values for primary proton, iron and 
mixed {p, He, O, Fe, Section 3.3) compositions computed in the framework of 
the SIBYLL and QGSJET interaction models. 



3 EAS inverse problem and primary energy spectra 

3.1 Key assumptions 

The observed spectra F{q) of the measured EAS parameters q = {Nch, N^, s) 
result from convolutions of the energy spectra Ia{E) of primary nuclei {A = 
H, He, ... at least up to Ni) with the probability density distributions Wa{E, q) 



The functions WA{E,q) are derived using a model of the EAS development 
in the atmosphere and convolution of the resulting shower spectra at the ob- 
servation level with the corresponding response functions [6,25]. 
The integral equation (1) defines the EAS inverse problem, namely the evalu- 
ation of the primary energy spectra Ia{E) on the basis of the measured distri- 
butions F{(\i) (in i = 1, . . . ,V discrete bins) and the known kernel functions 
WA{E,(\i) [6,25,40]. The multidimensional kernel functions Wa{E,c\) can be 
computed using interpolations [13] or approximations [6] to the corresponding 
spectra, which are previously obtained by CORSIKA EAS simulations in the 
framework of a given interaction model, for different groups of primary nuclei 
and a set of primary energies and zenith angles. 

In the present work, the computations of the expected shower spectra (forward 
folding) from (1) for given primary energy spectra Ia{E) are performed by 
Monte Carlo integration [42,43], using an arbitrary positive weight function 
Iq{A,E) determined in the same domain as the primary spectra Ia{E) and 
normalized such that / Wa{E)Iq{A, E) dE = 1. 

Multiplying and dividing the integrand in (1) by Io{A,E), expression (1) is 
converted to the form: 



The averaging in (2) is performed over random Ej (j = 1, . . . ,N'a) distributed 
with a probability density function Io{A,E)Wa{E \ qj e qj), with shower 



[13,33,40]: 




(1) 




(2) 
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parameters qj within the given bin. The reconstructed shower parame- 
ters q^j{A,Ej) are obtained by EAS simulations in the framework of a given 
interaction model, taking into account the corresponding response functions 
< SnJA^N,!,) > (Section 2.3). 

As a weight function we chose the power law spectrum Io{A, E) oc E~^-^ 
which provides an accuracy for integration AF/F ~ l/vCV" and relatively 
small statistical errors for the simulated EAS samples both within and espe- 
cially beyond the knee region. The accuracy of Monte-Carlo integration with 
this weight function was checked using power-law spectra f{x) oc x~'^ with 
7 = 2.5 — 3.3 and log-normal distributions W{x, y), and found to be adequate 
for our purposes. 

In order to evaluate the primary energy spectra on the basis of the EAS data 
set we regularized the integral equation (1) using a parametrization method 
[13,15]. The solutions for the primary energy spectra in (1) were sought based 

on a broken power-law function with a "knee" at the rigidity-dependent energy 
Ej:[A) = Er ■ Z, and the same spectral indices for all species of primary nuclei 
(^4 = p. He, O, Fe), 71 below and 72 above the knee respectively: 



where 7 = 71 for £" < Ek{A), 7 = 72 for £" > Ek{A), Er is the particle's 
magnetic rigidity and Z the charge of nucleus A. 

The integral equation (1) is thereby transformed into a parametric equation 
with the unknown spectral parameters $a, Er, 71 and 72, which are evaluated 
by minimization of the function: 



where \J is the number of examined functions Cm,^ = Fuisiu.i) obtained from the 
experimental data with statistical accuracies (j{C,u,i) in i = 1, . . . , bins, and 
^u,i and cr(^„,i) are the corresponding expected (forward folded) values from 
(2) and their (statistical) uncertainties. 



3.2 Simulated EAS database 



EAS simulations for the evaluation of the primary energy spectra using the 
GAMMA facility EAS data were carried out for A/a = 10^ primary H, 7.1 ■ 10"^ 
He, 4.6 • 10^ O and 4.8 • 10^ Fe nuclei using the CORSIKA NKG mode and 
the SIBYLL interaction model. The corresponding statistics for the QGSJET 
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interaction model were: 10^, 6 • 10^, 4.4 • 10^ and 4 • 10^. 

The energy thresholds of the primary nuclei were the same for both interaction 
models and were set at EA,rain = 0.5, 0.7, 1 and 1.2 PeV for H, He, O and Fe 
respectively, and the upper energy limit was set at i^max = 5-10^ PeV. The sim- 
ulated energies were distributed following a weight function Io{A, E) oc E~^-^, 
as explained above. The simulated showers had core coordinates distributed 
uniformly within a radius < 25 m, and zenith angles 9 < 30°. This ignores 
the effect of showers with true core coordinates outside the selection radius 
which have reconstructed coordinates with i? < 25 m; due to the good core 
reconstruction accuracy of 0.5 — 1 m (Section 2.2), this effect may be neglected 
for our purposes. 

All the EAS muons with energies of i?^ > 4 GeV at the GAMMA observation 
level were passed through the 2.3 kg/cm^ of rock to the muon scintillation 
carpet (the underground muon hall). The fluctuations in the muon ioniza- 
tion losses, and the electron (and positron) accompaniment due to the muon 
electromagnetic and photonuclear interactions in the rock were taken into ac- 
count, using the approximation of an equilibrium accompanying charged par- 
ticle spectrum obtained from preliminary simulations with the FLUKA code 
[41] in the 0.005 — 20 TeV muon energy range. The resulting charged particle 
accompaniment per EAS muon in the underground hall is equal to 0.06 ±0.01 
(100%e) and 11.0 ± 1.5 (98.5%e, 1.4%/i, 0.04%//) at muon energies 0.01 TeV 
and 10 TeV respectively. 

The total number of simulated EAS in the database were Af = I] A/a — 
2.65 ■ 10^ EAS for the SIBYLL and Af ~ 2.44 ■ 10^ EAS for the QGSJET 
model. 



3.3 Combined approximations to the EAS data 

Using the aforementioned formalism (Section 3.1), the U — 6 examined func- 
tions shown in Figs. 6-11 and the corresponding EAS data set, the unknown 
spectral parameters ^A,Ep>,^i and 72 of parametrization (3) were derived by 
minimization of the (4) and forward folding (2), with a number of degrees 
of freedom n^.f. — Y^tVu —p — 1 ~ 350, where p = 7 is the number of adjustable 
parameters. 

The values of the spectral parameters (3) derived from the solution of the 
parameterized equation (1) are presented in Table 1 for the SIBYLL and 
QGSJET interaction models. The primary energy spectra obtained for p. He, O, 
and Fe nuclei, along with the all-particle energy spectra, are shown in Fig. 12 
(hues and shaded areas) for the SIBYLL (left panel) and QGSJET (right 
panel) interaction models. The symbols in Fig. 12 show the all-particle spec- 
tra obtained by KASCADE [6] from a 2-dimensional {Ne, N^) unfolding using 
an iterative method, and from GAMMA [27] using an event-by-event method. 
Also shown as error bars in the left panel of Fig. 12 are extrapolations of 
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Table 1 

Parameters of the primary energy spectra (3) from combined approximations to 
the EAS data. The scale factors and particle rigidity respectively have units 
of (m2. s • sr ■ TeV)-^ and TV. 



Parameters 


SIBYLL 


QGSJET 


^ 12 


0.095 ± 0.008 


0.165 ± 0.005 


^ TT 


n inn + 019 


n non _|_ n nno 


$o 


0.034 ± 0.007 


0.008 ± 0.004 




0.024 ± 0.004 


0.013 ±0.005 


Er 


2500 lb 200 


3200 ± 150 


71 


2.68 ± 0.015 


2.66 ± 0.010 


72 


3.19 ± 0.03 


3.11 ±0.02 




2.0 


1.5 




Primary energy E (GeV) Primary energy E (GeV) 

Fig. 12. Energy spectra and abundances of the primary nuclei groups (lines and 
shaded areas) for the SIBYLL (left panel) and QGSJET (right panel) interaction 
models. All-particle spectra from GAMMA [27] and KASCADE [6] are shown as 
symbols. Vertical bars show the extrapolations of balloon and satellite data [37]. 



the balloon and satellite data to the energy E ~ 10^ GeV, computed us- 
ing power-law approximations to the available direct measurement data [37]; 
these remain in reasonable agreement with more recent balloon experiment 
data [38,39]. In this extrapolation, the 0-Iike group was assumed to include 
the elements Z — 3-16, and the Fe-hke group the elements Z — 17-28. 
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Tabic 2 

Parameters of the primary energy spectra (3) from 2-dimensional approximations 
to the EAS data. The scale factors and particle rigidity respectively have 
units of (m^- s • sr • TeV)~^ and TV. 



Parameters 


SIBYLL 


QGSJET 


^ 12 


0.109 ± 0.006 


0.198 ± 0.006 


^ TT 


f) nqc: _|_ n r\r\a 




$o 


0.050 ± 0.006 


0.031 ± 0.002 




0.017 ±0.002 


0.006 ± 0.002 


Er 


2500 lb 200 


4200 ± 300 


71 


2.70 ± 0.005 


2.71 ± 0.030 


72 


3.23 ± 0.08 


3.23 ± 0.09 


X^/nd.f. 


1.2 


1.3 



The expected EAS spectra and Nch{s) and NchiNj^) dependencies according 
to the solutions presented above are shown in Figs. 6-11 for the QGSJET 
(dashed hues) and SIBYLL (sohd lines and shaded areas) interaction models. 
The vertical widths of the shaded areas correspond to the error bars of the 
expected spectra, which are comparable for the two interaction models. 
It should be noted that the results obtained in the framework of the QGSJET 
interaction model strongly depend on the number of examined functions, 
which is not the case with the SIBYLL model. 



3.4 2- Dimensional approach 



Using (1), parametrization (3) and the 2-dimensional EAS spectra 



we evaluated the parameters of the primary energy spectra by minimization of 
the corresponding function (4), with U = 1. The computations were carried 
out with bin dimensions A\n Nch = 0.15 and AlnA*"^ = 0.25, for 6 < 30° and 
Nch > 5 ■ 10^. The resulting number of degrees of freedom {n^.f) for the 
minimization was equal to about 240. 

The best-fit spectral parameters and corresponding values of x^/n^./. for both 
interaction models are presented in Table 2. The contributions to the total 
from each 2-dimensional bin = {N^-h, N^) at the minimum of the x^ 
function are shown in Figs. 13 and 14, for the SIBYLL and QGSJET models 
respectively. 
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_ SIBYLL2.1 

7.5 



7 



6.5 



6 



3 3.5 4 4.5 5 

Lg(N^) 

Fig. 13. Contributions to the total from each ij^chi^ix) bin, for the SIBYLL 
interaction model. 



QGSJET01 




3 3.5 4 4.5 5 

Lg(Nti) 



Fig. 14. Same as Fig. 13 for the QGSJET interaction model. 
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Table 3 

Parameters of the primary energy spectra (3) from the 4-dimcnsional analysis of 
the EAS data. The scale factors and particle rigidity respectively have units 
of (m^- s • sr • TeV)-^ and TV. 



Parameters 


SIBYLL 


QGSJET 


^ 12 


0.110 ± 0.004 


0.190 ± 0.002 


^ TT 


f) nqi _|_ ri r\(\A 


n noo _|_ n nno 


$o 


0.045 ± 0.004 


0.038 ± 0.002 




0.030 ±0.002 


0.010 ± 0.002 


Er 


2300 lb 230 


3100 ± 200 


71 


2.67 ± 0.005 


2.68 ± 0.005 


72 


3.13 ± 0.06 


3.10 ±0.06 




2.1 


2.1 



3.5 4-Dimensional approach 



The amount of information about the primary energy spectra contained in the 
4-dimensional spectrum of measured parameters 



F(q) = 



dTVcft dNi^ ds d cos 9 



is obviously always greater than the information contained in the 2-dimensional 
{Nch, Nfj^) spectrum (Section 3.4) or the cumulative amount of information 
contained in the combined spectra (Section 3.3). The main difference with the 
latter case is due to the inter-correlations between EAS parameters, which can 
only be fully taken into account in such a 4-dimensional approach. 
On the basis of this 4-dimensional representation of the EAS data set, the 
simulated EAS database, and parameterization (3), equation (1) was solved 
by minimization, with U = 1. The computations were carried out with the 
following bin dimensions: Aln A^cft — 0.15, AlnA^^ = 0.25, Asec^ = 0.05, and 
As = 0.15 on the left and right hand side of s* — 0.85 and As = 0.3 else- 
where. The number of degrees of freedom in this 4-dimensional approximation 
was equal to 1640. The values of spectral parameters (3) resulting from the 
solution of the parameterized equation (1) are presented in Table 3 for the 
QGSJET and SIBYLL interaction models. 
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4 Discussion 



As can be seen from Fig. 12 and Tables 1-3, the derived primary energy spec- 
tra depend significantly on the interaction model, and slightly on the approach 
(Sections 3.2-3.5) applied to solve the EAS inverse problem. The derived abun- 
dances of primary nuclei at an energy E ~ 10^ TeV in the framework of the 
SIBYLL model agree (in the range of 1-2 standard errors) with the corre- 
sponding extrapolations of the balloon and satellite data [37], whereas the 
results derived with the QGSJET model point toward a dominantly proton 
primary composition in the 10^ — 10^ TeV energy range. 

Although the derived formal accuracies of the spectral parameters in Tables 1- 
3 are high, the corresponding values are large, which demands further dis- 
cussion. These large values do not necessarily imply disagreement between 
the EAS data and the derived primary energy spectra, but could be due to 
a number of other possible causes. We believe that the most likely causes of 
the large values of our spectral fits are systematic uncertainties related to 
the EAS simulations, in the interaction model or in the computation of the 
detector response (Section 2.3), and to the representation of the full cosmic 
ray composition by a small number of simulated nuclear species. We note that 
the inclusion of additional errors of about 5 — 7% in the functions (4) will 
decrease the x^/'^rf./. — 2 in Tables 1-3 to /ud.f. — 1. 

We discuss in turn below a number of other possible causes and related is- 
sues, especially the possibility that our spectral parametrization is incorrect, 
in terms of the rigidity-dependent knee energy or common spectral index. We 
also consider briefly the uncertainties in the reconstructed spectral parame- 
ters, discuss possible issues with the convergence of the unfolding method and 
the number of elemental groups, and present some consistency checks on the 
simulated and experimental databases. 



^.1 Rigidity-dependent knee hypothesis 

A test of the spectral parametrization (3) was performed by evaluating the 
knee energies Ek{A) independently for each primary nucleus A = H, He, O, Fe 
simultaneously with the spectral parameters ^AjTi and 72, using the combined 
approximation method described above (Section 3.3). The derived scale fac- 
tors and spectral indices 71 and 72 agreed with the data from Table 1 
within errors, but they had somewhat larger uncertainties (typically by fac- 
tors ~ 1.2 — 1.7). The derived knee energies versus nuclear charge (Z) are 
shown in Fig. 15 for the SIBYLL and QGSJET interaction models (symbols), 
along with the corresponding expected values (lines) according to the rigidity- 
dependent knee hypothesis from Table 1. It can be seen from Fig. 15 that the 
independently adjusted knee energies agree with the rigidity-dependent knee 
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Fig. 15. Knee energies for each group of nuclei (symbols) versus nuclear charge 

Z. The expected values in the rigidity-dependent (thick lines) and mass-dependent 
(thin lines) knee hypotheses are also shown, for the QGSJET (dashed lines) and 
SIBYLL (solid lines) interaction models. 



hypothesis. 

We also examined the alternative mass-dependent knee hypothesis; also shown 
as lines in Fig. 15 are the results of spectral fits using combined approxima- 
tions (Section 3.3), in which the hypothesis Ek{A) = Ek ■ A, with A the 
nuclear mass, was assumed. The values of the xLin were practically the same 
as in Table 1, but the derived value of the spectral parameter 71 tended to the 
range 2.59 ± 0.02, which is somewhat hard relative to expectations from the 
balloon and satellite data [37,38,39]. Within the uncertainties of our present 
analysis, our data are not in contradiction with this A-dependent knee hypoth- 
esis; however, it clearly does not yield a better agreement than our assumed 
rigidity-dependent hypothesis. 

4-2 Common spectral index hypothesis 

We attempted to similarly examine the possibility of independent spectral 
indices 71^^ for each primary nucleus, A = H, He, 0,Fe, but in that case en- 
counter a difficulty. The solution found by minimization when these param- 
eters are independent strongly depends on the initial values given to the min- 
imization algorithm, making a thorough exploration of the multi-dimensional 
parameter space impractical, and the results inconclusive. 
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Fig. 16. Xmin dependences versus the values of the fixed parameter (71) with ran- 
domly chosen initial values for the adjustable parameters (see text for details). 



Figure 16 shows the X^ili) dependences for different spectral hypotheses. The 
thick solid line represents the x^ili) for parameterization (3), where the spec- 
tral index is common to all primary species, obtained in the combined approx- 
imation approach (Section 3.3) with the SIBYLL interaction model. The other 
lines show the corresponding dependences x^{1i,a) individual nuclei, in the 
case where the lower spectral indices 71^^ are independent for each species. In 
all cases, the value of the parameter shown is held fixed, but values for all 
the other parameters are obtained by minimization of (4), with initial values 
for the minimization algorithm assigned randomly in a range of ~ 10 — 20% 
around representative values for the spectral parameters "Ji^a, 72 and logE'^, 
and in a range of ~ 50 — 100% for the scale factors It is readily seen that 
while the curve x^ili) for a common 71 shows a quite robust behavior, the 
minima found for independent spectral indices strongly depend on the initial 
values. 

The shape of the x^ili) curve for the parametrization with equal spectral 
indices (3) may be used as an illustration to examine the reliability of the un- 
certainties quoted in Tables 1-3. The minimization in all cases was performed 
using the FUMILI algorithm [45], and the errors quoted were obtained from 
the formal covariance matrix of the fit at the minimum. A more accurate es- 
timation of the confidence interval can be obtained from the intersection of the 
appropriate level Ax^ above the minimum value with a curve such as the 
thick solid line in Fig. 16. After normalizing the errors such that x^/n^./. ~ 1, 
we find that the actual confidence interval is somewhat wider than that ob- 
tained from the formal uncertainty. In general, our investigations suggest that 
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the derived formal errors tend to underestimate the actual uncertainties in the 
spectral parameters by up to a factor of two. 



4-3 Problem of uniqueness 



The example of independent spectral indices 71,^1 illustrates a more general 
potential difficulty. The EAS inverse problem is an ill-conditioned problem by 
definition, and unfolding of the corresponding integral equations (1) does not 
ensure the uniqueness of the solutions. Furthermore the EAS inverse problem 
implies the evaluation of two or more unknown primary energy spectra from 
an integral equation set of the Fredholm kind, and this peculiarity has not 
been studied in detail. 

Evidently, the solution cannot be considered unique if a small change in the 
initial values of the iterative algorithm for the minimization of (4) results 
in a significant change (well beyond the formal uncertainties) of the solution 
spectra. Using this test of uniqueness we concluded that only the equality 
of the spectral indices for all primary nuclei below the knee and the same 
equality of the spectral indices above the knee (parameterization (3)) result 
in the unique solutions presented in Fig. 12 and Tables 1-3. 



4-4 Number of elemental groups 



The evaluations of primary spectra for pure H, pure He and mixed {H, He), 
{H, He, O) and {H, He, Fe) compositions in parameterization (3) also were 
examined using the 2-dimensional approach (Section 3.4). The correspond- 
ing x^/'^d.f. values were respectively equal to 44.5,35.3,10.0,1.8 and 2.5 for 
the SIBYLL interaction model and 11.5, 141, 4.0, 2.7 and 2.0 for the QGSJET 
model. The results for mixed H, He, O and Fe primary composition are pre- 
sented in Table 2. It is readily seen that the data cannot be adequately rep- 
resented with less than the four considered types of primary nuclei. 
Examining these results we can conclude that increasing the number of con- 
sidered primary nuclei in our parameterized inverse approach increases the 
accuracy of the solutions. This effect indirectly supports the validity of our 
parametrization with equal spectral indices. If our assumption of the equality 
of the spectral indices was invalid, we would not expect the to improve so 
effectively with increasing number of nuclear species. 
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Fig. 17. EAS size spectrum with an enlarged shower core selection criterion 

(i? < 50 m) (symbols), and expected shower spectra for each of the primary nuclei 
and the mixed composition (lines and shaded areas), with parameters from Table 1 
and for the SIBYLL interaction model. 



4-5 Consistency of the solutions 



The agreement of the data presented in Tables 1 and 3 with our preliminary 
results [25,26], which were obtained with significantly fewer (half as many) 
simulated showers, suggests that the size of the simulated database is not a 
problem. 

A further check of the consistency of the GAMMA facility EAS data with the 
derived solutions is shown in Figs. 17 and 18, which display the EAS size and 
truncated muon size spectra (symbols) for an enlarged core selection criterion 
of i? < 50 m. This is twice as large as the selection radius of the EAS data 
in Figs. 6-7, and resulted in about four times the number of selected showers. 
The lines and shaded areas in Figs. 17-18 correspond to the expected (forward 
folded) EAS spectra with the parameters of primary energy spectra (3) from 
Table 1 for the SIBYLL interaction model; the corresponding expected shower 
spectra for each of primary nuclei are also shown. 
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Log,o(N^ (R^,<50m) ) 
Fig. 18. Same as Fig. 16 for the EAS truncated muon size spectra. 



5 Conclusion 

The consistency of the results obtained by the GAMMA experiment (Figs. 6- 
11,16-17), at least up to Ncu — 10^, with the corresponding predictions in 
the framework of the hypothesis of a rigidity-dependent knee in the primary 
energy spectra and the validity of the SIBYLL or QGSJET interaction models 
points toward the following conclusions: 

• A rigidity-dependent steepening of primary energy spectra in the knee re- 
gion (expression 3) describes the EAS data of the GAMMA experiment 
at least up to Nch — 10^ with an average accuracy < 10%, with particle 
magnetic rigidities Er ~ 2500 ± 200 TV (SIBYLL) and Er ~ 3100 - 4200 
TV (QGSJET). The corresponding spectral power-law indices are 71 = 
2.68 ± 0.02 and 72 = 3.10 — 3.23 below and above the knee respectively, and 
the element group scale factors are given in Tables 1-3. 

• The abundances and energy spectra obtained for primary p, He, 0-likc and 
Fe-like nuclei depend on the interaction model. The SIBYLL interaction 
model is preferable in terms of consistency of the extrapolations of the 
derived primary spectra (Fig. 12) with direct measurements in the energy 
range of satellite and balloon experiments [37,38,39] . 

• The derived all-particle energy spectra depend only weakly on the interac- 
tion model. They are compatible with independent measurements of this 
spectrum. 

• An anomalous behavior of the EAS muon size and density spectra (Fig. 5b, 
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Figs. 11 and 18) and the EAS age parameter (Fig. 10) for EAS size Nch > 10^ 
is observed. A similar behavior of the EAS age parameter has previously 
been observed in [30,44]. The observed behavior of the muon size and density 
spectra may be related to the excess of high- multiplicity cosmic muon events 
detected by the ALEPH and DELPHI experiments [46,47]. 
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